In [1]:
from __future__ import print_function, division
In [2]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
import time
import shutil
import urllib
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [3]:
import hosts
import decals
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord, Angle
from astropy import table
from astropy.table import Table
from astropy.io import fits
from astropy.utils import data
import tqdm
from IPython import display
In [4]:
%matplotlib inline
from matplotlib import style, pyplot as plt
plt.style.use('seaborn-deep')
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['figure.figsize'] = (14, 8)
In [5]:
bricks = Table.read('decals_dr4/survey-bricks.fits.gz')
bricksdr3 = Table.read('decals_dr3/survey-bricks-dr3.fits.gz')
bricksdr4 = Table.read('decals_dr4/survey-bricks-dr4.fits.gz')
In [6]:
paper1nsas_comp = [166313,147100,165536,61945,132339,149781,33446,150887]
paper1nsas_incomp = [161174,85746,145729,140594,126115,13927,137625,129237]
paper1nsas = paper1nsas_comp + paper1nsas_incomp
In [7]:
hostobjs = hosts.get_saga_hosts_from_google()
for host in hostobjs:
host.fnsdss = 'catalogs/base_sql_nsa{}.fits.gz'.format(host.nsaid)
hostsbyname = {h.name:h for h in hostobjs}
paperhosts = [h for h in hostobjs if h.nsaid in paper1nsas]
assert len(paperhosts) == len(paper1nsas)
In [8]:
host_bricks_3 = decals.find_host_bricks(paperhosts, bricksdr3, bricks)
host_bricks_4 = decals.find_host_bricks(paperhosts, bricksdr4, bricks)
host_bricks_dct = {3:host_bricks_3, 4:host_bricks_4}
In [9]:
for hostobj in paperhosts:
hbricks = host_bricks_3[host_bricks_3['closest_host_name'] == hostobj.name]
catfn = 'decals_dr3/'
In [10]:
def load_host_catalog(hostobj, drnum):
host_bricks = host_bricks_dct[drnum]
tabs = []
for brick in host_bricks[host_bricks['closest_host_name'] == hostobj.name]:
catfn = 'decals_dr{}/catalogs/tractor-{}.fits'.format(drnum, brick['brickname'])
cat = Table.read(catfn)
cat['objname'] = ['{}-{}'.format(row['brickname'], row['objid']) for row in cat]
tabs.append(cat)
if tabs:
return table.vstack(tabs)
In [11]:
for hostobj in tqdm.tqdm_notebook(paperhosts):
hostobj.cats3 = load_host_catalog(hostobj, 3)
[(hostobj.name, len(hostobj.cats3)) for hostobj in paperhosts if hostobj.cats3 is not None]
Out[11]:
In [12]:
for hostobj in tqdm.tqdm_notebook(paperhosts):
hostobj.cats4 = load_host_catalog(hostobj, 4)
[(hostobj.name, len(hostobj.cats4)) for hostobj in paperhosts if hostobj.cats4 is not None]
Out[12]:
In [13]:
for host in tqdm.tqdm_notebook(paperhosts):
for cat in (host.cats3, host.cats4):
if cat is not None:
decals.mags_catalog(cat)
decals.aperture_sbs_catalog(cat, bandname='r')
decals.interpolate_catalog_sb(cat, loopfunc=lambda x: tqdm.tqdm_notebook(x, leave=False))
In [14]:
psfdepths = []
galdepths = []
for host in paperhosts:
for cat in (host.cats3, host.cats4):
if cat is not None:
scs = SkyCoord(cat['ra'], cat['dec'], unit=u.deg)
subcat = cat[host.within_environs(scs)]
if 'decam_galdepth' in subcat.colnames:
psfdepthr = subcat['decam_depth'][:, 2]
galdepthr = subcat['decam_galdepth'][:, 2]
else:
psfdepthr = subcat['psfdepth_r']
galdepthr = subcat['galdepth_r']
psfdepths.append(-2.5*(np.log10(5*psfdepthr**-0.5)-9))
galdepths.append(-2.5*(np.log10(5*galdepthr**-0.5)-9))
In [15]:
psfdepthsall = np.concatenate(psfdepths)
galdepthsall = np.concatenate(galdepths)
plt.hist(psfdepthsall[np.isfinite(psfdepthsall)], bins=100, histtype='step', label='psf')
plt.hist(galdepthsall[np.isfinite(galdepthsall)], bins=100, histtype='step', label='gal')
plt.legend(loc=0)
None
In [16]:
host = hostobjs[0]
cat = hostobjs[0].cats3
scs = SkyCoord(cat['ra'], cat['dec'], unit=u.deg)
sep = scs.separation(host.coords)
cat = cat[sep < host.environsarcmin*u.arcmin]
In [17]:
fig, ax1 = plt.subplots(1,1)
psfs = cat['type'] == 'PSF '
lowsb_saga = (cat['decam_mag'][:, 2] < 21) & (cat['sbeff_r']>24.5) & np.isfinite(cat['sbeff_r']) & ~psfs
print('nlowsb=', np.sum(lowsb_saga))
ax1.scatter(cat['mag_r'][psfs], cat['sbeff_r'][psfs], alpha=.15, lw=0, c='r')
ax1.scatter(cat['mag_r'][~psfs], cat['sbeff_r'][~psfs], alpha=.25, lw=0)
ax1.scatter(cat['mag_r'][lowsb_saga], cat['sbeff_r'][lowsb_saga], alpha=.35, lw=0, c='g')
ax1.set_ylabel('sberr')
ax1.set_xlim(15, 26)
ax1.set_ylim(15, 30)
decals.show_decals_objects_in_nb(cat[lowsb_saga], 4, info_cols=['mag_r', 'sbeff_r'])
Out[17]:
In [18]:
fig, ax1 = plt.subplots(1,1)
cat['decam_mag_r'] = cat['decam_mag'][:, 2]
psfs = cat['type'] == 'PSF '
lowsb_saga = (cat['decam_mag'][:, 2] < 21) & (cat['sb_r_0p5']>24.5) & np.isfinite(cat['sb_r_0p5']) & ~psfs
print('nlowsb=', np.sum(lowsb_saga))
ax1.scatter(cat['mag_r'][psfs], cat['sb_r_0p5'][psfs], alpha=.15, lw=0, c='r')
ax1.scatter(cat['mag_r'][~psfs], cat['sb_r_0p5'][~psfs], alpha=.25, lw=0)
ax1.scatter(cat['mag_r'][lowsb_saga], cat['sb_r_0p5'][lowsb_saga], alpha=.35, lw=0, c='g')
ax1.set_ylabel('sberr')
ax1.set_xlim(15, 26)
ax1.set_ylim(15, 30)
decals.show_decals_objects_in_nb(cat[lowsb_saga], 4, info_cols=['decam_mag_r', 'sb_r_0p5'])
Out[18]:
This strongly suggests that sbeff_r is the way to go
In [19]:
for host in tqdm.tqdm_notebook(paperhosts):
try:
sdsscat = host.get_sdss_catalog()
except IOError:
print('Missing base catalog for', host.name)
sdsscat = None
for cat in (host.cats3, host.cats4):
if cat is not None:
d2ds = np.full(len(cat), -1, dtype=float)
if sdsscat is not None:
dcatsc = SkyCoord(cat['ra'], cat['dec'], unit=u.deg)
idx, d2d, _ = dcatsc.match_to_catalog_sky(sdsscat['coord'])
plt.figure()
plt.hist(d2d.arcsec, bins=1000, histtype='step',range=(0, 30), log=True)
plt.axvline(1, c='k', ls='--')
plt.axvline(3, c='k')
plt.title(host.name)
d2ds = d2d.to(u.arcsec)
cat['dist_basecat'] = d2ds
cat['dist_basecat'].description = '-1 means "unknown", e.g. base catalog missing'
In [ ]:
In [51]:
maxnamelen = max([len(h.name) for h in paperhosts])
tostack = []
for host in tqdm.tqdm_notebook(paperhosts):
for drnum, cat in [(3, host.cats3), (4, host.cats4)]:
if cat is not None:
psfs = cat['type'] == 'PSF '
lowsb_saga = (cat['mag_r'] < 21) & (cat['sbeff_r']>24.5) & np.isfinite(cat['sbeff_r']) & ~psfs
subcat = cat[lowsb_saga]['objname', 'ra', 'dec', 'mag_r', 'sbeff_r', 'dist_basecat']
subcat['dr'] = np.full(len(subcat), drnum, dtype=int)
subcat['hostname'] = np.full(len(subcat), host.name, dtype='S'+str(maxnamelen))
tostack.append(subcat)
towrite = table.vstack(tostack)
towrite.sort(['hostname', 'dist_basecat'])
towrite.meta = {}
towrite['inbasecat'] = np.full(len(towrite), 'unknown', dtype='S'+str(len('unknown')))
towrite['inbasecat'][towrite['dist_basecat'] > 0*u.arcsec] = 'yes'
towrite['inbasecat'][towrite['dist_basecat'] > 1*u.arcsec] = 'maybe'
towrite['inbasecat'][towrite['dist_basecat'] > 3*u.arcsec] = 'no'
towrite['inbasecat'].description = 'estimate of "is this in the base catalogs": <1 arcsec is yes, 1-3 is maybe, >3 is no'
towrite
Out[51]:
In [58]:
htmlyes = decals.show_decals_objects_in_nb(towrite[towrite['inbasecat']=='yes'], dr='fromcatalog', nrows=4,
info_cols=['mag_r', 'sbeff_r', 'hostname', 'dr', 'dist_basecat'], sdss_link=True)
htmlnotyes = decals.show_decals_objects_in_nb(towrite[towrite['inbasecat']!='yes'], dr='fromcatalog', nrows=4,
info_cols=['mag_r', 'sbeff_r', 'hostname', 'dr', 'dist_basecat'], sdss_link=True)
with open('lowsb_decals_notinbasecat.html', 'w') as f:
f.write('<!doctype html>\n<meta charset=utf-8>\n<title>LowSB SAGA Decals object not in basecat</title>\n')
f.write('<h1>Non-basecat overlaps</h1>\n')
f.write(htmlnotyes.data)
with open('lowsb_decals_inbasecat.html', 'w') as f:
f.write('<!doctype html>\n<meta charset=utf-8>\n<title>LowSB SAGA Decals objects in basecat</title>\n')
f.write('<h1>Already in basecats</h1>\n')
f.write(htmlyes.data)
In [59]:
towrite.write('/Users/erik/Dropbox/SAGA/temporary/lowsbs_decals.dat', format='ascii.ecsv')
!cp lowsb_decals*.html /Users/erik/Dropbox/SAGA/temporary/